{
 "cells": [
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T05:53:33.958382Z",
     "start_time": "2025-07-31T05:53:33.687212Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import lstsq\n",
    "path = \"D:/7_30data/\""
   ],
   "id": "7641de60b7e629ba",
   "outputs": [],
   "execution_count": 2
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T07:04:00.613677Z",
     "start_time": "2025-07-31T07:04:00.603543Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def iterative_refinement_solve(V, L_rho, max_iter=25, tolerance=1e-15, verbose=True):\n",
    "    \"\"\"\n",
    "    高精度迭代精化求解复数线性系统 S @ a = b，其中 S = V^H @ V, b = V^H @ L_rho\n",
    "    专门追求机器精度，允许更多迭代次数，更严格的收敛标准\n",
    "\n",
    "    Parameters:\n",
    "    -----------\n",
    "    V : ndarray, shape (m, n)\n",
    "        输入复数矩阵\n",
    "    L_rho : ndarray, shape (m,)\n",
    "        右端复数向量\n",
    "    max_iter : int, default=25\n",
    "        最大迭代次数（增加到25次以确保充分迭代）\n",
    "    tolerance : float, default=1e-15\n",
    "        收敛容差（机器精度级别）\n",
    "    verbose : bool, default=True\n",
    "        是否输出详细信息\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    a : ndarray\n",
    "        解向量（如果虚部可忽略则返回实数，否则返回复数）\n",
    "    \"\"\"\n",
    "\n",
    "    # 确保输入数据类型为complex128以获得最高精度\n",
    "    V = V.astype(np.complex128)\n",
    "    L_rho = L_rho.astype(np.complex128)\n",
    "\n",
    "    # 构建正规方程 S @ a = b\n",
    "    S = V.conj().T @ V  # V^H @ V\n",
    "    b = V.conj().T @ L_rho  # V^H @ L_rho\n",
    "\n",
    "    # 计算矩阵的秩\n",
    "    rank_S = np.linalg.matrix_rank(S)\n",
    "    augmented = np.hstack((S, b.reshape(-1, 1)))\n",
    "    rank_augmented = np.linalg.matrix_rank(augmented)\n",
    "    print(f\"rank(S) = {rank_S}\")\n",
    "    print(f\"rank([S | b]) = {rank_augmented}\")\n",
    "\n",
    "    if verbose:\n",
    "        # 输出条件数信息用于诊断\n",
    "        cond_S = np.linalg.cond(S)\n",
    "        print(f\"Condition number of S: {cond_S:.4e}\")\n",
    "        print(\"\\n=== High Precision Iterative Refinement ===\")\n",
    "\n",
    "    # 检查系统的相容性\n",
    "    if rank_S != rank_augmented:\n",
    "        if verbose:\n",
    "            print(\"Warning: System is inconsistent (rank(S) ≠ rank([S|b]))\")\n",
    "            print(\"Using least squares solution...\")\n",
    "\n",
    "    # 步骤1: 获得初始解（使用最高精度的最小二乘法）\n",
    "    try:\n",
    "        a = lstsq(S, b, lapack_driver='gelsy')[0]\n",
    "        if verbose:\n",
    "            print(\"Initial solution computed using high-precision least squares\")\n",
    "    except Exception as e:\n",
    "        if verbose:\n",
    "            print(f\"Failed to compute initial solution: {e}\")\n",
    "        return np.zeros(S.shape[1], dtype=np.complex128)\n",
    "\n",
    "    # 计算初始残差\n",
    "    initial_residual = np.linalg.norm(S @ a - b)\n",
    "    if verbose:\n",
    "        print(f\"Initial residual norm: {initial_residual:.16e}\")\n",
    "\n",
    "    # 迭代精化过程 - 追求机器精度，不做过早的停滞检测\n",
    "    best_residual = initial_residual\n",
    "    best_solution = a.copy()\n",
    "    consecutive_no_improvement = 0\n",
    "\n",
    "    for i in range(max_iter):\n",
    "        # 步骤2: 计算当前残差 r = b - S @ a\n",
    "        r = b - S @ a\n",
    "        residual_norm = np.linalg.norm(r)\n",
    "\n",
    "        if verbose:\n",
    "            print(f\"Iteration {i+1:2d}: Residual norm = {residual_norm:.16e}\")\n",
    "\n",
    "        # 记录最佳解\n",
    "        if residual_norm < best_residual:\n",
    "            best_residual = residual_norm\n",
    "            best_solution = a.copy()\n",
    "            consecutive_no_improvement = 0\n",
    "        else:\n",
    "            consecutive_no_improvement += 1\n",
    "\n",
    "        # 步骤3: 严格的收敛检查\n",
    "        if residual_norm <= tolerance:\n",
    "            if verbose:\n",
    "                print(f\"Converged to tolerance after {i+1} iterations\")\n",
    "            break\n",
    "\n",
    "        # 步骤4: 求解修正方程 S @ da = r\n",
    "        try:\n",
    "            da = lstsq(S, r, lapack_driver='gelsy')[0]\n",
    "        except Exception as e:\n",
    "            if verbose:\n",
    "                print(f\"Iteration {i+1} failed: {e}\")\n",
    "            break\n",
    "\n",
    "        # 步骤5: 更新解 a = a + da\n",
    "        a = a + da\n",
    "\n",
    "        # 只在连续多次无改善且远未达到机器精度时才考虑停止\n",
    "        # 这里设置更严格的条件：连续8次无改善且残差仍大于10倍机器精度\n",
    "        if (consecutive_no_improvement >= 8 and\n",
    "            best_residual > 10 * tolerance and\n",
    "            i >= 10):  # 至少迭代10次\n",
    "            if verbose:\n",
    "                print(f\"Stopping: No improvement for {consecutive_no_improvement} consecutive iterations\")\n",
    "                print(\"Using best solution found so far\")\n",
    "            a = best_solution.copy()\n",
    "            break\n",
    "\n",
    "    # 计算最终残差（使用最佳解）\n",
    "    final_residual = np.linalg.norm(S @ a - b)\n",
    "\n",
    "    if verbose:\n",
    "        print(f\"\\nFinal residual norm for S @ a - b: {final_residual:.16e}\")\n",
    "\n",
    "        # 输出改进情况\n",
    "        improvement_factor = initial_residual / final_residual if final_residual > 0 else float('inf')\n",
    "        print(f\"Improvement factor: {improvement_factor:.2e}\")\n",
    "\n",
    "        # 检查解的性质\n",
    "        max_imag = np.max(np.abs(np.imag(a)))\n",
    "        if max_imag > 1e-12:\n",
    "            print(f\"Solution has complex components (max imaginary part: {max_imag:.4e})\")\n",
    "        else:\n",
    "            print(\"Solution is effectively real (converting to real)\")\n",
    "\n",
    "    # 输出必要的最终信息（不截断，显示完整精度）\n",
    "    print(f\"Residual norm for S @ a - b: {final_residual:.16e}\")\n",
    "\n",
    "    # 智能返回类型：如果虚部可忽略则返回实数\n",
    "    max_imag = np.max(np.abs(np.imag(a)))\n",
    "    if max_imag < 1e-12:\n",
    "        return np.real(a)\n",
    "    else:\n",
    "        return a"
   ],
   "id": "1834cae8b1f232c3",
   "outputs": [],
   "execution_count": 34
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T07:04:10.572981Z",
     "start_time": "2025-07-31T07:04:09.527144Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 使用示例_1\n",
    "V = np.load(path+'V.npy').astype(np.complex128)\n",
    "L_rho = np.load(path+'L_rho.npy').astype(np.complex128)\n",
    "x = iterative_refinement_solve(V, L_rho,verbose=True)\n",
    "print(\"\\n\")\n",
    "y = iterative_refinement_solve(V, L_rho,verbose=False)"
   ],
   "id": "7922e930f19b4780",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Condition number of S: 9.7122e+80\n",
      "\n",
      "=== High Precision Iterative Refinement ===\n",
      "Initial solution computed using high-precision least squares\n",
      "Initial residual norm: 1.4461906702066530e-10\n",
      "Iteration  1: Residual norm = 1.4461906702066530e-10\n",
      "Iteration  2: Residual norm = 2.0579515875677572e-11\n",
      "Iteration  3: Residual norm = 8.2318063498103511e-11\n",
      "Iteration  4: Residual norm = 2.0579515875517370e-11\n",
      "Iteration  5: Residual norm = 2.0579515875476124e-11\n",
      "Iteration  6: Residual norm = 2.0579515875463642e-11\n",
      "Iteration  7: Residual norm = 2.1084465186937517e-16\n",
      "Converged to tolerance after 7 iterations\n",
      "\n",
      "Final residual norm for S @ a - b: 2.1084465186937517e-16\n",
      "Improvement factor: 6.86e+05\n",
      "Solution is effectively real (converting to real)\n",
      "Residual norm for S @ a - b: 2.1084465186937517e-16\n",
      "\n",
      "\n",
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Residual norm for S @ a - b: 2.1084465186937517e-16\n"
     ]
    }
   ],
   "execution_count": 35
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T06:57:14.489111Z",
     "start_time": "2025-07-31T06:57:11.272593Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 使用示例_2\n",
    "np.random.seed(42)\n",
    "V = np.random.randn(4, 273) + 1j * np.random.randn(4, 273)\n",
    "L_rho = np.random.randn(4) + 1j * np.random.randn(4)\n",
    "solution = iterative_refinement_solve(V, L_rho, verbose=True)\n",
    "solution = iterative_refinement_solve(V, L_rho, verbose=False)\n",
    "\n",
    "# 使用示例_3\n",
    "np.random.seed(100)\n",
    "V = np.random.randn(4, 273) + 1j * np.random.randn(4, 273)\n",
    "L_rho = np.random.randn(4) + 1j * np.random.randn(4)\n",
    "solution = iterative_refinement_solve(V, L_rho, verbose=True)\n",
    "solution = iterative_refinement_solve(V, L_rho, verbose=False)"
   ],
   "id": "c4bde89c257229a9",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Condition number of S: 1.8539e+19\n",
      "\n",
      "=== High Precision Iterative Refinement ===\n",
      "Initial solution computed using high-precision least squares\n",
      "Initial residual norm: 2.8866405161597366e-14\n",
      "Iteration  1: Residual norm = 2.8866405161597366e-14\n",
      "Iteration  2: Residual norm = 1.7600757088117352e-14\n",
      "Iteration  3: Residual norm = 1.7711692175792960e-14\n",
      "Iteration  4: Residual norm = 1.7279235285330439e-14\n",
      "Iteration  5: Residual norm = 1.8574849046989312e-14\n",
      "Iteration  6: Residual norm = 1.9666570987670721e-14\n",
      "Iteration  7: Residual norm = 1.8971336933154305e-14\n",
      "Iteration  8: Residual norm = 1.8019034268696099e-14\n",
      "Iteration  9: Residual norm = 1.7507053463922958e-14\n",
      "Iteration 10: Residual norm = 1.8728776492714919e-14\n",
      "Iteration 11: Residual norm = 1.8043750063518024e-14\n",
      "Iteration 12: Residual norm = 1.7860327410391217e-14\n",
      "Stopping: No improvement for 8 consecutive iterations\n",
      "Using best solution found so far\n",
      "\n",
      "Final residual norm for S @ a - b: 1.7279235285330439e-14\n",
      "Improvement factor: 1.67e+00\n",
      "Solution has complex components (max imaginary part: 3.1180e-02)\n",
      "Residual norm for S @ a - b: 1.7279235285330439e-14\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.7279235285330439e-14\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Condition number of S: 1.0462e+19\n",
      "\n",
      "=== High Precision Iterative Refinement ===\n",
      "Initial solution computed using high-precision least squares\n",
      "Initial residual norm: 4.5320149444698993e-14\n",
      "Iteration  1: Residual norm = 4.5320149444698993e-14\n",
      "Iteration  2: Residual norm = 1.4480890553958694e-14\n",
      "Iteration  3: Residual norm = 1.6180686109998086e-14\n",
      "Iteration  4: Residual norm = 1.5660338464658046e-14\n",
      "Iteration  5: Residual norm = 1.6769108501221586e-14\n",
      "Iteration  6: Residual norm = 1.5028491273321592e-14\n",
      "Iteration  7: Residual norm = 1.5817913333999005e-14\n",
      "Iteration  8: Residual norm = 1.5588640973570810e-14\n",
      "Iteration  9: Residual norm = 1.5697245042087597e-14\n",
      "Iteration 10: Residual norm = 1.5809097817073499e-14\n",
      "Iteration 11: Residual norm = 1.5670046012356502e-14\n",
      "Stopping: No improvement for 9 consecutive iterations\n",
      "Using best solution found so far\n",
      "\n",
      "Final residual norm for S @ a - b: 1.4480890553958694e-14\n",
      "Improvement factor: 3.13e+00\n",
      "Solution has complex components (max imaginary part: 3.8662e-02)\n",
      "Residual norm for S @ a - b: 1.4480890553958694e-14\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.4480890553958694e-14\n"
     ]
    }
   ],
   "execution_count": 32
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T06:58:52.368039Z",
     "start_time": "2025-07-31T06:58:26.875217Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 测试\n",
    "# 设置测试组数量\n",
    "num_test_cases = 20\n",
    "\n",
    "# 使用循环一次生成多组数据并测试函数\n",
    "for i in range(num_test_cases):\n",
    "    print(f\"--- test {i+1} ---\")\n",
    "    # 为每个测试用例设置不同的随机数种子，以确保数据独立且可复现\n",
    "    rng = np.random.default_rng(seed=57 + i)\n",
    "    # 使用当前的 rng 对象生成 V 和 L_rho\n",
    "    V = rng.standard_normal(size=(4, 273))+ 1j * rng.standard_normal(size=(4, 273))\n",
    "    L_rho = rng.standard_normal(size=(4)) + 1j * rng.standard_normal(size=(4))\n",
    "    # 调用函数并测试\n",
    "    try:\n",
    "        solution_silent = iterative_refinement_solve(V, L_rho, verbose=False)\n",
    "    except Exception as e:\n",
    "        print(f\"  test {i+1} ，error : {e}\")\n",
    "    print(\"-\" * 30 + \"\\n\")"
   ],
   "id": "2d3538714533e34b",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- test 1 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.2743260353250655e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 2 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 8.7977275403987738e-15\n",
      "------------------------------\n",
      "\n",
      "--- test 3 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.7417833211758851e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 4 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.0744090478420666e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 5 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.9385244789896381e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 6 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.0067997907628746e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 7 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.2900318220698719e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 8 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.0090702013698455e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 9 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.0217485341986212e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 10 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.1630243025309659e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 11 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.3759902937319692e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 12 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.5267787365929146e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 13 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.1040891333171151e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 14 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.3609700699767815e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 15 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.2186699476798661e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 16 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.1390105572420087e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 17 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.2038496896953706e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 18 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.2346784876108291e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 19 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.4484541409240567e-14\n",
      "------------------------------\n",
      "\n",
      "--- test 20 ---\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.0207704804375586e-14\n",
      "------------------------------\n",
      "\n"
     ]
    }
   ],
   "execution_count": 33
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}
